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Abstract. We describe the computational techniques employed in the re- 
cently updated Fourier local correlation tracking (FLCT) method. The FLCT 
code is then evaluated using a series of simple, 2D, known flow patterns that 
test its accuracy and characterize its errors. 



1. Introduction - What is the basic concept behind local correlation 
tracking (LCT)? 

Given a pair of 2D maps ("images") h(x, y, t) and h(x, y, t2 = tl + St) of some 
scalar quantity, with the second image taken slightly later than the first one, 
what is the 2-d flow field (y x [x, y], v y [x, y]) which, when applied to the scalar 
field in the first image, will most closely resemble the second image? This 
definition of LCT is not a precise one, and the LCT technique incorporates no 
physical conservation laws. Schuck (2005) showed that LCT methods are more 
consistent with an advection equation, rather than a continuity equation. Here, 
we acknowledge these shortcomings and forge ahead, noting that LCT results 
must be interpreted carefully. 

In Solar Physics, the idea for LCT is generally attributed to November and 
Simon (1988). In the engineering literature, the problem is known as the "optical 
flow" problem (see Schuck [2006] and references therein). Here, we present 
the technique behind the recently upgraded Fourier LCT (FLCT) method, first 
described in Welsch, Fisher, & Abbett (2004). 

2. The Mathematical Approach Used by FLCT 

To construct a 2D velocity field that connects two images Ii(x,y) and l2(x,y) 
taken at two different times t\ and ti-, one must start from some given location 
within both images, compute a velocity vector, and then repeat the calculation 
while varying that location over all pixel positions. This involves three high-level 
operations: (1) windowing the input images to isolate the neighborhood around 
the pixel of interest; (2) computing the correlation function between the two 
images; and (3) locating the peak of the cross correlation function. 

For each pixel at which a velocity is to be computed, a windowing function 
is used to de-emphasize parts of the image far away from that pixel. FLCT 
does this localization by multiplying each of the two images by a Gaussian of 
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width a, centered at pixel location (x{,yj). We denote the resulting images as 
"sub-images" Si and S2. The expressions for Si and S2 are: 

S$ J) (x,y) = I 2 {x,y)e^ x - x ^ + (y-y^ 2 l° 2 . (1) 

The quantity a is a free parameter in FLCT, and its optimal value changes 
depending on the nature of the image and the size scales present in the velocity 
field. As discussed in Welsch et al. (2007), one way to choose the optimal a for 
a given application is to select the a that statistically minimizes the difference 
between the final image and the advected initial image, \Ii — 6t(y ■ V)ii|. 

For the (i,j)th pixel, the cross-correlation function of sub-image 1 with 
sub-image 2 is defined by 

C^(6x,Sy) = J J dxdySi' j *(-x,-y)S i /(5x-x,5y-y) . (2) 

We want to find, for each pair of sub-images Si and S2 centered at position 
(xi,yj), the shifts 5x and 5y that maximize C(5x,5y). The amplitude of the 
shifts, divided by the time 5t = ti — t\ between images 1 and 2 defines the 
velocity determined by FLCT: v x = 5x/5t, and v y = Sy/St. 

FLCT uses the convolution theorem to compute C(5x,5y) using Fourier 
transforms. If we write J 7 (Si) = si(k x ,k y ) and, T(S2) = S2(k x ,k y ), where T 
denotes Fourier transform, then the above equation can be written 

C i j{5x,5y)=T- l {s\s 2 ) , (3) 

where J-^ 1 denotes the inverse Fourier transform. We sometimes find it useful 
to perform a low-pass Gaussian filter on the functions s± and S2 before applying 
equation @ if the original images are noisy. Other researchers have used differ- 
ent cross-correlation functions; some calculate the correlation in x,y-space, e.g., 
November and Simon (1988). 

Next, we must find the peak of the cross-correlation function. Sub-pixel 
accuracy is required, as shifts are frequently substantially less than 0.1 pixel. 
As a practical matter, we find the peak of f(Sx,8y) = \C l ' :, {5x, Sy)\, rather 
than C(5x,5y) so that the operation does not involve complex arithmetic. For 
notational simplicity, we henceforth use x, y for <5x, 5y in the followin g discussion. 



Previous versions of FLCT followed ( November and Simon 1988T ). and inter- 



polated / around its pixel-accuracy peak onto a finer grid, and took the location 
of the maximum of the resulting discretely sampled function as the shift. This 
approach, however, is only as accurate as the resolution of the interpolated grid, 
and is therefore computationally expensive — many unnecessary interpolations 
are required to reach the necessary spatial resolution. 

Our current version employs a curve-fitting approach to find the peak in 
f(x,y) that was inspired by that of Chae (2004, private communication; LCT 
code written in IDL). First, since the images and sub-images are computed at 
discrete points in space, we identify the pixel coordinates (x m , y n ) of the largest 
value of /, denoting the largest value of / as f(x m ,y n ). Note that (x m ,y n ) m ay 
not be equal to (xi, yj), if the location of the peak of f(x, y) has shifted by more 
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than a pixel in x or y. To find the peak to sub-pixel resolution, we then Taylor- 
expand f(x, y) to 2nd order about the (x m ,y n ) location, denoting the expansion 
as fr(x,y), 

Of df 

fr{x,y) = f(x m ,y n ) + — (x - x m ) + —{y - y n ) (4) 

ox oy 

+ ^y {x ~ Xm){v ~ Vn) + 2^ {X ~ Xm) + 2W {V ~ Vn) ' 

where the partial derivatives are evaluated at the point (x m ,y n ). 

At the peak, we require that the x and y partial derivatives of the Taylor 
expansion vanish. These conditions result in a pair of linear equations 

which allow us to solve for the location (x max ,y max ) of the peak: 



-i 

(5) 



(d 2 fdf d 2 f df\ (d 2 f\ 2 d*fd*f 

\dy 2 dx dxdy dy J \\dxdy I dx 2 dy 2 

(d 2 fdf d 2 f df\((d 2 f\ 2 d 2 fd 2 f\~ 1 

I dx 2 dy dxdy dx J \ \dxdy I dx 2 dy 2 J 

To evaluate the partial derivatives, we use standard, second-order finite differ- 
ence expressions, assuming a uniform grid in the sub- images: 

df df 

-q^ = (/m+l,n — /m-l,n)/2Al , — = (f m ,n+l ~ /m,n-l)/2Ay 

d 2 f 2 ® 2 f 2 

— (/m+l,n frn—l,n ^frn,n)/^X , Qy2 — {fm,n+l fm,n— 1 2/ mn )/Ay 

d 2 f 

r, = (/m+l,n+l + fm-l,n-l ~ fm-l,n+l ~ /m+l,n-l)/4AxAl/ . (7) 

dxdy 

The total shift in (x, y) that corresponds to the peak in / is then given by 

5x = (x m - Xi) + (x max - Xm) , 5y = (y n - yj) + (y max - y n ) . (8) 

Equations (|5|) through ([8]) thus provide a simple and accurate method for finding 
the peak of the cross- correlation function to sub-pixel resolution. 



3. The Computational Approach Used by FLCT 

Our method requires two input images, I\ and I2, recorded at times t± and 
£2 = ti + dt. The algorithm then loops over all pixels in 1\ and I2 for which 
\h +I2 1/2 > hhr-, where I t h r is a user-set threshold. For each such pixel, it must: 

1. Compute sub- images S\ and S2. 

2. Compute Fourier transforms of S\ and S2, si and 82- 

3. Perform low-pass filters on si, S2 if needed. 
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4. Compute inverse Fourier transform of s*«2- This is C i,J . 

5. Compute the absolute value of C l > 3 , \C l ' J \. 

6. Compute the shifts 5x,8y that maximize |C l,:, |. 

7. Compute velocities v x = 5x/5t,v y = Sy/6t. 

The FLCT code, which is currently called flct, was initially written in 
IDL, but has since been re-written in C for portability and speed. FLCT also 
uses its own endian-independent binary input format for the images and for the 
output velocity arrays. FLCT uses the FFTW library (version 3) for computing 
the Fast-Fourier Transforms needed to compute the cross-correlation functions. 
We have written IDL procedures which write the two images flct needs into 
an input file, and which read the flct output file, consisting of two-dimensional, 
floating point arrays of v x ,v y , and a mask array, v m , equal to one in pixels with 
\Ii + -^2 1/2 > ithr and zero elsewhere. These IDL i/o utilities enable flct to be 
run easily from within an IDL session. 

4. Running FLCT 

To compile the FLCT code, the only external library needed is the FFTW3 
library. To download a copy of the FLCT source code and compilation instruc- 
tions, along with IDL i/o procedures, go to our web siteQ Be sure to get the C 
version (not the IDL version), and get version test_13 or later. 

The flct is code invoked from the command- line - either a terminal window 
on a Unix-like machine, from the command-prompt tool in MS Windows, or 
from within an IDL session using either a shell escape character or the spawn 
command. In all cases, the syntax for flct is as follows: 

flct infile outfile deltat deltas sigma -t thr -k kr -h -q 

Required Arguments: 

infile - Contains 2 images for local correlation tracking. To create infile 
use the IDL procedure vcimage2out .pro 

outfile - contains output v x ,v y ,v m (x, y velocity and mask arrays). Mask 
array is zero where velocity not computed. To read outfile, use the IDL 
procedure vcimage3in.pro 

deltat - Amount of time between images. 

deltas - Units of length of the side of a single pixel; velocity is computed 
in units of deltas/deltat. 

sigma - Sub-images are weighted by Gaussian of width sigma. If sigma 
is set to 0, only single values of shifts are returned as v x ,v y . These values 
correspond to the overall shifts between the two images, and these values 
can be used to co-register the images if desired. 
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Optional parameters: 

-t thr - Determines threshold parameter I t hr- If \h + ^2 1/2 < hhx in a 
pixel, then skip calculation of shifts for that pixel. For thr € (0, 1), thr 
is assumed to be in relative units of the maximum absolute pixel value in 
{|/i|, 1^2 1 }• To force thr values between zero and unity to be in absolute 
units, append an 'a' to thr. If calculation is skipped for a pixel, then the 
v m mask array is set to zero at that pixel location. It is otherwise set to 
1. 

-k kr - Apply a low-pass filter to the sub- images, with a Gaussian of a 
characteristic wavenumbcr that is a factor of k r times the largest possible 
wave numbers in x,y directions. k r should be positive. (If kr >> 1, the 
unfiltered result should be recovered). This option is sometimes useful for 
datasets that contain noise at high spatial frequencies. Test cases in which 
MDI Quiet-Sun magnetograms are shifted by 10 milli-pixels frequently 
result in a factor of two under-estimate by f let of the applied shift. By 
setting k r to numbers between 0.2 to 0.5, the full applied shifts can be 
recovered. Experimentation is recommended. 

-h - "hires" mode. Set this flag to use cubic convolution interpolation of 
the f(x, y) function to re-grid it to higher resolution before finding shift. It 
is unclear that this interpolation improves accuracy. Currently, this option 
seldom used, and may be dropped in the future. It has been retained for 
heritage. 

-q - Flag to suppress printing of all non-error messages. 

From within IDL, a simple demonstration of FLCT can be run using the 
following commands. 

IDL> fl=randomu(seed,101,101) 

IDL> f2=shift(f 1,1,-1) 

IDL> vcimage2out , f 1 , f 2 , ' testin . dat ' 

IDL> $flct testin.dat testout.dat 1. 1. 15. 

IDL> vcimage3in , vx , vy , vm , ' testout . dat ' 

IDL> shade_surf , vx 

IDL> shade_surf , vy 



5. FLCT Speed 

How quickly does FLCT run? On a 1.6 GHz Pentium M (laptop processor) 
running MS Windows XP, the FLCT code Act can process roughly 200 pixels 
per second for images with over 40,000 pixels (time per pixel is much less than 
this for smaller images). The code is speeded up dramatically if the -t thr 
option is used, in which the velocity calculation is skipped if the image values 
are less than the chosen threshold. There are many avenues available for further 
speedup of the FLCT code. Very little experimentation has been done with 
multiple threads thus far. An experimental version of Act written in Fortran, 
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using MPI has been tested on our Linux cluster, grizzly, and the compute speed 
was found to scale almost perfectly with the number of processors. 




Figure 1. An MDI magnetogram of the Quiet Sun. By applying known 
shifts and rotations to such magnetograms, we can test the performance of 
FLCT. 



6. Tests of FLCT 



Figure Q] shows an MDI ( Scherrer et al. 19951 ) magnetogram of the Quiet Sun. 



(For a discussion about applying tracking techniques like FLCT to magne- 
tograms, see the paper Welsch and Fisher in this volume.) To test the ability of 
FLCT to reconstruct rotations, we rotated the image, using cubic-convolution 
interpolation, by a fixed angle about the center of the image to generate a sec- 
ond image. The applied "velocity field" then corresponds to simple, solid-body 
rotation, with the speed increasing linearly with radius, r, from image center. 

How well does FLCT recover the applied velocity field? In Figures [2] and O 
we show scatter plots of derived speeds versus input speeds for a rotations of 1° 
and 0.01°, respectively. Units of velocity are in pixel widths per unit time. The 
straight lines show the applied speed, and the scatter plots show the inverted 
speeds for all inverted pixels (only points with \B\ > 10 G were inverted), a = 
15 pixels was assumed. For the 1° rotation, the scatter is typically 0.1-0.2, 
and there is no significant bias in the derived speed. For the 0.01° rotation, 
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Figure 2. A scatter plot of derived speeds versus input speeds, for an ap- 
plied rotation of 1°. The straight line shows the applied speed. There is no 
significant bias in the derived speed. 



the scatter is significant, and the derived speeds are biased on the low side as 
compared to the applied speeds. For the 1° rotation, the recovered rotation 
profile (not shown) accurately reproduces the imposed rotation profile. These 
and other tests suggest that FLCT can accurately recover shifts > 0.1 pixels, 
but smaller shifts can be problematic. 

7. Discussion 

What is the difference between the new (upgraded) FLCT code, described here, 
and that originally described in Welsch et al. (2004)? The main difference is 
in how the location of the peak of the cross-correlation function is determined 
to sub-pixel resolution. The old version used cubic convolution interpolation to 
compute the cross-correlation function on a much finer grid (0.1 or 0.02 pixel 
spacing) and then simply found the location of the largest value within the finer 
grid. We (and other users) found that this resulted in serious "quantization" 
errors when the time between images was small enough that shifts were of order 
0.1-0.2 pixels or less. The new code uses equations ([7]) and (jSJ) of this paper 
to find the location of the peak to sub-pixel resolution. We have found this to 
be more accurate, and more computationally efficient — and therefore much 
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Figure 3. A scatter plot of derived speeds versus input speeds, for an applied 
rotation of 0.01°. The scatter is significant, and the derived speeds are biased 
on the low side as compared to the applied speeds. 



faster — than the old technique. We strongly recommend that users use the 
new version and not any of the old versions. 
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